Fluctuation-regularized Front Propagation Dynamics 
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We introduce and study a new class of fronts in finite particle number reaction-diffusion systems, 
corresponding to propagating up a reaction rate gradient. We show that these systems have no 
traditional mean-field limit, as the nature of the long-time front solution in the stochastic process 
differs essentially from that obtained by solving the mean-field deterministic reaction-diffusion equa- 
tions. Instead, one can incorporate some aspects of the fluctuations via introducing a density cutoff. 
Using this method, we derive analytic expressions for the front velocity dependence on bulk particle 
density and show self-consistently why this cutoff approach can get the correct leading-order physics. 



Many physical, chemical, and biological systems ex- 
hibit fronts which propagate through space. Familiar 
examples range from chemical reaction dynamics such 
as flame fronts [ij, phase transitions such as solidifica- 
tion 2] , the spatial spread of infections 3] , and even the 
fixation of a beneficial allele in a population Ijj. In all 
these cases, the underlying dynamics of individual con- 
stituents (from molecules to organisms) gives rise mul- 
tiple macroscopic states. It is thus of great interest to 
understand the universality classes of fronts which gov- 
ern what will happen when systems such as these are pre- 
pared in a spatially heterogeneous manner. These classes 
determine the selection of propagation speed, the sensi- 
tivity to particle-number fluctuations, and the stability 
of the front with respect to deviations from planarity. 

To date, several different classes of fronts have been 
described. Perhaps the simplest one, exemplified by the 
case of solidification, is that wherein a thermodynami- 
cally stable phase replaces a metastable one 2]. Here 
the mean-field front velocity is determined via the re- 
quirement that there exists a heteroclinic trajectory of 
the steady-state problem (obtained by assuming that all 
fields depend only on the traveling coordinate z = x — vt) 
connecting the metastable phase at +cx3 with the stable 
one at — oo. This type of front is robust with respect to 
fiuctuations, with power-law corrections in 1/A^ (where 
A'' is the number of particles per site in the final state) 
to the mean- field limit [a|- A second class is exemplified 
by the simple infection model A + B ^ 2A on Id lattice 
with equal A and B hopping rates 3f|; this process leads 
in the mean-field limit to the well-known Fisher equa- 
tion y (^t = r0(l — (j)) + D(j)x:x- Here propagation is 
into the linearly unstable = state and (p is just the 
number of A particles at a site, normalized by A'^. Re- 
cent work [a 0> Si has shown that the front behavior 
in the stochastic model does approach that of the Fisher 
equation, where the velocity is selected by the (linear) 
marginal stability criterion [2j to be 2-\/rD, albeit with 
an anomalously long transient 0{l/t) and anomalously 



large fiuctuation corrections 0(1/ In A^). There are also 
some surprising findings in regard to both front stability 
in the case of unequal D llfl j, a nd also the scaling prop- 
erties of front fiuctuations [ll|. Finally, there are also 
fronts which are determined by the nonlinear marginal 
stability criterion (for example, propagation into a non- 
linearly unstable but linearly marginal state) which have 
properties intermediate to the previous two classes. 

In this work, we introduce a new class of fronts corre- 
sponding to propagation up a reaction-rate gradient. We 
focus again on the A + B ^i- 2 A reaction '3| , with no A 
particles and an initial mean number A^ of B particles at 
all sites past some initial xo, but now assume that the 
reaction probability depends on spatial position. The 
two situations we wish to consider are respectively the 
absolute gradient and the quasistatic gradient 



raix) = max(rm„,ro-Hax) , 
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instantaneous position 
we identify as Xf = 



where Xf{t) is the 
of the front, which 

lJ2i=o^Aixi)xi/J2t=oNA{xt), where on the lat- 
tice Xi = ia. (The minimum value of r is just there to 
ensure that the bulk remains stable for all x and plays 
no important role.) The quasistatic problem should 
lead to a translation-invariant front with fixed speed 
Vq{'ro,a), whereas the absolute gradient will lead to an 
accelerating front. In the latter case, one can imagine 
ignoring the acceleration, obtaining thereby an adiabatic 
approximation to the velocity v{t]ro,a) ~ Vq{fo{t),a) 
where fo(i) — tq + axf{t). As we will see, fiuctuation 
effects are absolutely crucial as the naive mean-field 
theory gives rise to "irregular" fronts in a way which 
we will define shortly. It should be noted that glimpses 
of this new class were obtained in studies of models 
of Darwinian evolution [l^ |l3, |lj|, but no general 
understanding was attained. 

In Fig. 1, we show numerical solutions of the mean 
field equation (MFE, here just the Fisher equation with 
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FIG. 1: Comparison of the MFE equation and the stochastic 
model, for D = 5, ro = 1 and a — .02. (a) Numerical inte- 
gration of the MFE, with initial conditions </!> = 1 for x < xo, 
otherwise, (b) Twenty simulations of the stochastic system 
with A'^ = 10®, for each of two xq's, xq = 10, 40; also plot- 
ted is the corresponding cutoff MFE graph with k = 0.3. We 
note in passing that the distribution of velocities seems highly 
skewed, with a large tail at velocities in excess of the mean. 



spatially varying r) versus a stochastic stimulation |5lll5l| 
of the original continuous time Markov-process, both for 
the absolute gradient case. The MFE behaves irregularly 
in the sense that the front never recovers from its initial 
conditions, i.e. the system never falls into a dynamic 
attractor. Conversely, we define a regular front as one 
for which changes in initial data (as long as Na remains 
zero past some starting point) can only effect a time- 
translation of an otherwise fixed front solution and hence 
would be invisible on a plot of velocity versus position. 
Clearly, the stochastic process gives rise to a regular front 
and thus cannot in any manner be approximated by the 
MFE. 

To get some insight into the notion of regularity, we 
employ a heuristic approach in which we mimic the 
leading-order effect of finite population number fluctu- 
ations by introducing a cutoff in the MFE [13, [H, llTJ . 
This cutoff replaces r{x) by zero if the density (j) falls 
below k/N for some 0(1) constant k; this change in the 
reaction term prevents the leading edge from spreading 
too far, too fast. This idea has proven its reliability in the 
Fisher system with constant reaction rate where it cor- 
rectly predicts the aforementioned anomalous effects |6|. 
A numerical simulation of the cutoff equation reveals the 
recovery of regularity (see Fig. 2a) in its long-time be- 
havior; of course, the time it takes to converge to the 
dynamic front attractor diverges as A^ ^ oo (data not 
shown). Returning to the simulation results in Fig lb, 
we see that the cutoff MFE (using fc as a fitting param- 
eter) does a quantitatively accurate job of tracking the 
actual front dynamics. Fig. 2b presents a comparison 
of the front profile from the cutoff theory with that of 
the stochastic model; this was done for the quasistatic 
model as the translation invariance facilitates the neces- 
sary ensemble averaging. Overall, it is quite remarkable 
how well this simple approach works; later, we will use 
our analytic approach to the cutoff MFE to give at least 
some new indications of why this might be the case. 



FIG. 2: (a) Regularity of the cutoff MFE; same parameters 
as Fig. la except for the cutoff of 10^''. (b) Average profile 
from twenty simulations of the (quasistatic) stochastic system 
with A*' — lO'^; also plotted is the corresponding cutoff MFE 
graph with k — 0.3 



We next turn to using the cutoff MFE to study the 
front velocity. In Fig. 3, we present results for the 
front velocity at a fixed spatial position for the stochastic 
model and the cutoff mean field theory as a function of 
A^. This is done both for the absolute gradient model 
and for the corresponding quasistatic model. From the 
data, we can draw the following conclusions. First, both 
models exhibit velocities which increase without bound 
as a function of A^; this is of course radically different 
than what has been encountered in the previous classes 
of propagating fronts. This behavior accounts for the 
fact that the long-time dynamics never approaches that 
predicted by the naive mean-field theory. Next, the two 
different models exhibit similar velocities at small log A^ 
but become increasingly different as log A' goes to in- 
finity. Finally, we note that at small enough log A', the 
velocity can be approximated by just taking a cutoff ver- 
sion of the usual Fisher equation result for a fixed reaction 
rate rp = rg+ax, i.e. completely neglecting the reaction- 
rate gradient across the front. This can be explained by 
noting that the effective interfacial width, the distance 
over which the particle density drops from its bulk value 
0(1) to its cutoff value 0{1/N) scales as log A^; hence 
one can neglect the gradient if alog A^ is small. 

Given the highly unusual velocity results, an analytic 
treatment of the cutoff MFE at large A^ is clearly worth- 
while. We restrict our attention to the quasistatic case 
where the problem is a "standard" one of finding a homo- 
clinic trajectory in the traveling coordinate z. In what 
follows, it will become clear that spatial discretization 
effects cannot be neglected and hence we keep the ex- 
plicitly discrete form of the hopping term, with lattice 
spacing a. The first key idea is that the nonlinearity is 
important only near the bulk state and in that region, 
diffusion is much less important than drift if the velocity 
is large. In this range of sites, then, the full equation 



D((/)(z + a)+ (j){z -a)~ 2(f){z)) 
+v4>' + r(z)(0 - (f)^)0{4) - k/N) 
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FIG. 3: Cutoflt MFE equation for D = 5, ro = 1 and a = 0.1, 
velocity at a; = 200 versus logN. Shown for comparison is 
the absolute gradient model, the quasistatic model with fo = 
ro + 200q, a numerical solution of a cutoff version of the usual 
Fisher equation with r/ = ro + 200a and the Brunet-Derrida 
analytic approximation thereof. 



can be solved by neglecting diffusion; this leads to 
1 — (l)\ ( rQ{x + ax^/2)/v 



In 



rmin{x - x^)/v + rQ{x^ + ax^^/2)/v 



for X ^ Xni respectively, where Xm is the point where 
the minimum reaction rate is reached, i.e. rg + axm = 
Tni. We have chosen (^(0) = 1/2, to fix the translation 
invariance, which is a slightly different definition of the 
front location, amounting to a small shift in rg. 

The above result needs to be matched to the solu- 
tion near the cutoff point, where the nonlinearity can be 
dropped. As initially suggested by Rouzine, et al. [l4|. 
the resultant linear equation can be approached via the 
WKB ansatz (j) = e'^^^-', with the key observation that 
the cutoff point must occur close to the WKB turning 
point where (j) begins to oscillate. This is the only way 
that there is enough freedom to do the matching. The 
WKB equation takes the form 







AD 



sinh^(aS'72) + vS' + r^ + az 



(3) 



Already at this point, we get a nontrivial result. We 
can de-dimensionalize this equation by introducing T = 

a/iS, y = z/i, £ = v/{aa) so that the equation reads 



0= —smh^{r/2) + T' 
va 



y 



(4) 



where the derivative is now w.r.t. y and we have dropped 
the small term roa/ii.jl^l Thus, S (i.e. IniV, once we do 
the matching) scales like D/{aa?) times a function of the 
dimensionless parameter va/ D^ so that the leading-order 



results for all a, (for a given a and D) should lie on a uni- 
versal curve. Furthermore, we recover the idea already 
discussed in Ref. ([lj|) that a is a singular perturbation 
as far as the large velocity limit goes, since no matter 
how small a is, the parameter vajD eventually goes to 
infinity. 

Returning to Eq. Q , the turning point is given by the 
discriminant equation ^ = 0, yielding 



2D 



sinh(aS'^) -f v , 



which gives 



Sr = In 




2D 



(5) 



(6) 



If we match to the bulk solution near z = 0, (f> declines 
by an amount related to the change in S from z = to 
the turning point 2* . This is given by 



On: 



Si 



dz 



'''''-dS' 

S't f]Q' on 

— S'[-i—smh{aS')+v)] 



-^S'^ cosh(aS';|,) sinh(aS';|,) + -{Slf 



In the geometrical optics approximation, k/N « 0(zc) ~ 
0(z*) — e *, giving us an our lowest-order answer for 
the velocity. To improve upon this answer, one needs 
to both go beyond the geometrical optics approximation, 
and to develop a connection formula to go past the turn- 



ing point. The latter involves writing 



=Slx 



-0, where 



tp is smooth on the lattice scale [IJl, and showing that 
tfj satisfies an Airy equation. This procedure will be de- 
scribed in detail in a future publication 20], and here we 
just cite the results. The most significant correction, we 
find, comes from the further exponential decay, at rate 
SI, of (f> between the turning point and the cutoff point, 
which is near the zero of the Airy function solution to the 
ip equation in the connection region. This gives an ad- 
ditional contribution of —Sl^Q{k/ {D cosh{aSl)))~^/^ to 
S, where ^0 = —2.3381 is the location of the first zero 
of the Airy function. As already mentioned, the solution 
for the velocity is just ln(iV/fc) = —S, so that the final 
expression becomes 



ln(iV/fc) = - 
a 



a^ a-^ 2 



'S,oSn: 



Dcosh{aS'^] 



l3 

-1/3 



(7) 



Let us examine the various limits of the this expres- 
sion. First, in the continuum limit, av/D <ti 1. Then 
Si = v/2D, so that aS' < 1. Then, S*, = {2D{Sif/3 + 




MFE approach, we can understand in detail the anoma- 
lous velocity behavior, at least for the quasistatic case. 
Further work should address the role of acceleration and 
front stability. 

The work of HL has been supported in part by the 
NSF-sponsored Center for Theoretical Biological Physics 
(grant numbers PHY-0216576 and PHY-0225630). 



FIG. 4: (a) Velocity, v, vs. ln(A'^) for lattice spacings a = 1, 
0.5, 0.25, and 0, a = 0.1, D = ro = 1. from simulation, 
together with the analytic approximation Eq. (|7J (b) Velocity, 
u, vs. Q at In A'^ = 25, versus leading-order and full analytic 
expressions. In both figures, fc = 1. 



v{Si)'^/2)/a; the correction term is -Sl^o{k/D)-^^^. 
Combining these gives the formula 

V = 2(D2^oa)'/^ [(31n(iV/fc))i/3 +eo(31n(iV/fc))-i/3" . 

(8) 
The scaling of v as In^' '^ N is an agreement with the 
results of Tsimring, et al. for a continuum evolu- 
tion model [l3|. Now, as we mentioned above, for fi- 
nite a, av/D is eventually large for sufficiently large 
N. Then, Si = -\n{va/ D)/a. This gives S^ = 
v/{a^k){-\n{va/D) + 1 + \n{va/D)/2). Now, for very 
large va/D, S^ ~ —vh\{va/D)/{a?k). However, this is 
only valid for ln{va/D) ^2. In fact, it is a reasonable 
(20%) approximation only for \n{va/ D) bigger than 10, 
so that V would be unreasonably large. Thus, a strict 
asymptotic expansion is of no use. Our analytic results 
based on equationHare compared to numerical solutions 
in Fig. 4; more details are given elsewhere 20]. Note that 
the full expression fits the data quite well, although the 
leading-order formula is not very accurate. 

One of the most interesting aspects of the above cal- 
culation is that it does at all depend on form of the solu- 
tion past the cutoff point; the mere existence of a cutoff 
is enough to force the system to the WKB turning point 
and hence fix the velocity. This is perhaps the reason why 
the cutoff MFE mimics the actual stochastic model; the 
fact that the region past the cutoff is in reality highly 
stochastic should not alter the velocity fixation which 
occurs in the part of space where fluctuation effects are 
indeed negligible. Turning this qualitative argument into 
a full theory remains a challenge for future work. 

In summary, we have introduced a new class of fronts in 
reaction-diffusion systems and showed how fluctuations 
must be taken into account, even if only heuristically, 
if one wished understand their behavior. Using a cutoff 
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